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Abstract 

This note proposes rapidly convergent computational formu¬ 
lae for evaluating scattering kernels from radiative transfer 
theory. The approach used here does not rely on Legen¬ 
dre expansions, but rather uses exponentially convergent nu¬ 
merical integration rules. A closed form for the Henyey- 
Greenstein scattering kernel in terms of complete elliptic 
integrals is also derived. 
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1 Introduction and Background 

In theories of radiative transfer and of neutron transport, 
the interaction between radiation and a scattering medium is 
described by a phase function p : [—1,1] -A R + with the 
property p(t)dt = 1. Let S 2 denote the unit sphere 
in R 3 . In a single scattering event, a photon or neutron 
that arrives from a given direction © £ S 2 is scattered into 
the direction ©' £ S 2 according to the probability density 
©' i-A 4 ^rp(© • ©0- This probability density then appears 
as an integral kernel in an integro-differential equation that 
describes multiple scattering. In the special case of plane- 
parallel scattering, one rewrites this equation in spherical co¬ 
ordinates and expands it as a Fourier series in the azimuthal 
angle This requires the evaluation of 

n2n 

Cm / p(cos 6 > cos /i + sin 6 sin fi- cos (/)') cos me/) d(j)' ( 1 ) 

Jo 

for 0 < #, n < 7r, m — 0, 1,2,..., where c m = ^ for m > 0 
and co = Using the notation x = cos#, y = cos/q one 
therefore arrives at the problem of evaluating the scattering 
kernels P m (aj, y), defined for — 1 < x, y < 1 and for m — 
0, 1 ,... as 

p2tv 

Cm / p(xy + y/l — x 2 ^/l — y 2 cos s) cos ms ds (2) 

Jo 

Let L n denote the n-th Legendre polynomial. From spherical 
harmonics one obtains that the Legendre expansion of the 


phase function p directly leads to expansions for the P m . If 
p(x ) — a nL n (x ) for all x £ [—1,1], then in particular 

oo 

Po(x,y) — ^ ot n L n (x)L n (y) (3) 

71 = 0 

for all — 1 < x, y < 1. The higher order terms P m , m > 0 can 
be expanded into associated Legendre functions, using again 
the Legendre coefficients of p. The classical reference [2] con¬ 
tains a mathematical presentation of the theory of radiative 
transfer. A modern account with more physical details may 
be found in 6:. 

While the evaluation of eq. § and its higher order ver¬ 
sions is in theory straightforward, there are several practical 
difficulties. Firstly, the formula requires knowledge of the 
Legendre expansion of p. While the methods to find p for 
a particular scattering medium indeed produce Legendre ex¬ 
pansions (see [0), it may be desirable to have more compact 
representations of a scattering function, in which case the 
Legendre expansion is not readily available and not easy to 
compute (there is no “fast Legendre transform”). Secondly, 
for cases of strongly forward-peaked scattering, several hun¬ 
dred terms of the Legendre expansion may be needed to eval¬ 
uate each Pm even to modest accuracy. This requires care 
when evaluating the Legendre polynomials in eq. ([3]), and it 
is computationally expensive in any case. 

In this note, the direct numerical integration scheme known 
as trapezoidal rule is proposed to evaluate P m from eq. 
as an alternative to the Legendre expansion in eq. . The 
method is known (E3) to converge exponentially if the inte¬ 
grand is analytic. This highly desirable property is exploited 
systematically in this note. The relation between the con¬ 
vergence rate and the location of the singularities (points or 
regions of non-analyticity) of the phase function is explained. 
The second contribution of this note is the derivation of a 
closed form of the scattering kernel Po, in the special case of 
the Henyey-Greenstein phase function ([5]). It expresses Po 
in terms of a complete elliptic integral and can be evaluated 
very rapidly without any expansions or numerical integra¬ 
tion, even for cases of extremely forward-peaked scattering. 
Numerical examples are given to demonstrate the approach. 
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1.1 Henyey-Greenstein Scattering Func¬ 
tion 

The Henyey-Greenstein phase function (0)Q was proposed 
to describe interstellar scattering and is given by the formula 


PHa{x,g) 
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where — 1 < g < 1 is known as the asymmetry factor. This 
phase function has since been used in areas as diverse as scat¬ 
tering in cloudy and hazy atmospheres ( 0 ), light scattering 
in seawater ( 0 ) and in tissue ( 0 ), and even in computer 
graphics. Then eq. © together with eq. © lead to the prob¬ 
lem of evaluating 


°° o p I 1 

H(x,y;g) = ^2 —^— g l Li{x)L t {y) (5) 

1=0 

for —1 < x, y < 1. The series may be evaluated by using the 
first N terms. Since \Li{x)Li(y)\ ~ 21+1 indeterminate 
sign, the direct evaluation of the series leads to problems 
when g is close to 1, because then the series converges very 
slowly. This is illustrated in fig. [T] 


£71 



Figure 2: Log Errors for evaluating Henyey-Greenstein 
scattering kernel in fig. 0 Blue: Eq. (|5| with N = 40 
terms. Black, red, green: Eq. (19) with N = 40, 80, 160 
terms. Thick: Computed errors. Thin: Error estimates 


from eq. (20) 
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Figure 1: Henyey-Greenstein scattering kernel 

H(x,y 0 ;g) for g = 0.95, 2/0 = 0.4, -1 < x < 1. 
Black: exact evaluation using eq. (15). Blue: Eq. © 
with N = 40 terms. Red: Eq. (19) with N = 40 terms. 


1 It appears that Chandrasekhar was not aware of this work 
when he wrote his classical treatise [2 1 in 1950. 


2 An Exact Formula 


We start with the product formula (see 18.17.6 in m for 
Legendre polynomials 

1 r 

L n (cos 0)L n (cos n) = — L n (cos#cos/x 
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+ sin 6 sin g cos s)ds . 

Let 

00 

H 0 (x,y;g) = E g l L e (x)L e (y). (6) 

1=0 

Using the generating function for Legendre polynomials (see 
18.12.11 in m , we therefore obtain 


74o(cos 6 , cos g] g) 

1 PTV 00 

= — g n L n (cos# cos/i + sin 0 sin/x cos s) ds 
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ds 
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yjl — 2g (cos 6 cos g + sin 6 sin g cos s) + g 2 


(7) 

(8) 

(9) 


By eq. (33) with 
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a + [3 


1 + g 2 — 2g cos 0 cos g 
2g sin 6 sin g 
1 + g 2 — 2g cos (6 + g) 
























































this implies 


Hq (cos 6 , cos fi] g) = 


Ko 


7 + J3 + ^ 
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( 10 ) 


where Ao is the complete elliptic integral of the first kind 
defined in eq. (28). Setting 


u± = 1 + g 2 — 2g cos (0 ± /x) 
this may also be written as 


(ii) 


Ho (cos 0 , cos /x; g) = 
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XX+ — XX_ 
U+ 


This formula is closely related to formula (5.10.2. lQ in [9j 
which has the form 


Ho (cos 6 , cos /x; g) — 
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fv^+_ 


\Vu+ + V^ 


From eq. |l0l and eq. (p2| we obtain finally 


H (cos 0 , cos /x; g ) 
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where Eo is the complete elliptic integral of the second kind 
defined in eq. ([30|. In terms of x, y , this becomes 


H{x, y\g) = (1 p= E 0 

7 TW- ■ S /W+ 


4gVl - a; 2 a /1 - y 2 
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with w± = 1 + g 2 — 2g (xy =F \f\ — * 2 yi — y 2 ). The for- 
mula may be used to evaluate if(-, -;g) reliably even if g is 
extremely close to 1. If g = 1 — £ and £o ~ 10 -16 is “machine 
epsilon” in IEEE arithmetic, then H(x , y\g ) can be evaluated 
to relative accuracy £o/s. 


3 Fast Evaluation 

We now turn to the general case. For a given phase function 
p and given x, y G [— 1,1], m G {0, 1,2,...}, we need to 
evaluate the integral given by eq. © . Let 


the strip S a = {z G C | \5sz\ < a} and satisfies \f(z)\ < M 
for some constant M > 0 there. Choose a positive integer N 
and set 

' 2nk\ 
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then 
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f(t)dt-l N \ < eaJV _ 1 
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The sum is just an approximation of the integral with the 
composite trapezoidal rule. The rate of convergence thus is 
much faster than for ordinary smooth (twice differentiable) 
functions, for which the error estimate has the form 




" im = 


( 21 ) 


for some £ G (0, 2tt). 

To use this result in the computation of eq. ([5]), note first 
that for any qGR, the function that takes u + i • v = z to 

A + B cos z — A -\- B cos u cosh v + B sin u sinh v • i 


maps the strip S a to an ellipse about [A — B, A + B] in the 
complex plane which has focal points A ± B and major and 
minor axes with lengths Bcosha and B sinh a. Therefore, if 
the phase function p is analytic in a neighborhood surround¬ 
ing the set [—1,1] C C, then is analytic in a suitable strip 
S a with a > 0. The domain of analyticity of p is automat¬ 
ically symmetric with respect to the real axis. It should be 
emphasized that it is of course not necessary to determine 
this domain in order to use eq. |l9|). 

For an illustration, refer to fig. [3] The plot shows the do¬ 
main of an assumed phase function that is originally defined 
on the interval [—1,1] (thin horizontal black line) and that 
can be extended into the complex plane (everywhere except 
at singularities shown as colored lines and circles). For a 
particular choice of x , y, the integral in eq. © extends over 
[A — B, A + B] C [—1,1] (black circles). A suitable strip S a 
is mapped to the ellipse surrounding this set where the phase 
function is analytic. Consequently, the trapezoidal approxi¬ 
mation converges exponentially with a rate given by eq. (20). 


The rate of convergence depends on a which in turn comes 
from the location of [A — B,A + B] relative to the set of 
singularities of the phase function. 


hm(z) = p(A + B cos z) cos mz (17) 

A — xy , B — \J\ — x^yJX — y 1 . (18) 

Note that |Ad=5| < 1, with equality in one of the cases where 
x — ±y. The function hm is 27r-periodic on R. It is known 
(US) that if a function / is periodic and analytic in a strip 
about the real axis, then the trapezoidal rule for approxi¬ 
mating the integral J Q 27r f(t)dt converges exponentially fast. 
More precisely, assume that / is 27r-periodic and analytic in 

2 This was pointed out to me by an anonymous contributor on 
math.stackexchange.com. 


Example 1 Henyey-Greenstein Phase Function. Con¬ 
sider a phase function that is analytic in C minus the ray 
,oo). Then the integrand in eq. is analytic in the 
strip defined by 


cosh *Az < 


(i + 9 2 )/2g — A 


(l+g 2 )/2g -xy 

B 


Vl -x 2 ^l- y 2 


This case is illustrated by the magenta line in fig. [3] Similar 
comments apply to the case where g < 0 (red line). The 
Henyey-Greenstein phase function defined in eq. @ is of this 




































form. The error behavior of the trapezoidal rule ( |19| ) is illus¬ 
trated in fig. [5] which shows the log-errors (thick lines) when 
this rule is used to approximate the integral in eq. § for 
x G [— 1 , 1 ] and for fixed y — yo, g , for three different choices 
of N. The plot shows that errors are maximal for x ~ yo and 
that the errors have the slowest decrease near this value as 
N increases. Also plotted (thin lines, same colors) are error 
estimates from eq. (20), using M — 1 for simplicity. It is seen 
that the actual errors track the estimate very closely. 


Example 2 Consider a phase function that is analytic ex¬ 
cept possibly on rays xo ±i5 to xo d=z 00 , where xo E R, 5 > 0 . 
Then the integrand in eq. 0 is analytic in the strip defined 
by 

— S < sinh Sz<S ( 22 ) 

where S is the unique positive solution of the equation 


(A-xof f 2 
S 2 +1 S 2 


To see this, find a such that the ellipse parametrized by 

s 1 —»• A + B( cos s cosh a + i • sin s sinh a) 

passes through the points xq ± i • 5, and set S = sinh a. 

To obtain examples, fix S > 0, xo £ R, 7 > 0 and consider 
functions of x £ R 


where the proportionality constants are chosen such that the 
integrals over [—1,1] equal 1. Both functions have single 
maxima (peaks) at x = xo and the width of the peak is pro¬ 
portional to 5. The Legendre expansions of these functions 
are generally not available in closed form. The function fi is 
analytic in the complex plane minus branch cuts from xq d=z 5 
to xo =bz cx) (green lines in fig.[3|. The function f<i is analytic in 
the complex plane minus poles at z — xo±i8 (7t/2 + nir ), n £ 
Z + (blue circles in fig.p). Therefore if p = fi or p = / 2 , then 
the integrand in eq. ( 12 ) is analytic in any strip S a where 
a = arsinh S and S is as in example [ 2 ] 

The reader may note that the integrand in eq. contains 
factors cos mz. These terms grow like away from the 

real axis and their second derivatives contain the factor m 2 . 


When the integrand is analytic and eq. (20) can be used, the 
error is therefore proportional to e -aiV e am M = 

Thus the additional factor cos mz has the same effect as using 


N—m points instead of N points for the evaluation of eq. (19), 
resulting in a modest loss of accuracy. On the other hand, 
if t he in tegrand is merely twice differentiable, the error from 
eq. (21) becomes proportional to m 2 /N 2 . Thus the additional 
factor cos mz now has the same effect as using only N/m 
points instead of N points, leading to a much larger loss of 
accuracy. This illustrates the powerful effect of having an 
analytic integrand. 


3.1 Multimodal Phase Functions 


h(x\x 0 ,m,i) 

- (■*(?)!’ 

(23) 

f 2 (x;x 0 ,Tn) 

(X - Xo\ 
oc sech —-— J 

(24) 


In practice, phase functions are obtained from scattering cal¬ 
culations using Mie theory, see e.g. m Such phase functions 
may have multiple local extrema. An artificial example (not 
obtained from Mie theory) is given in fig. [4] It uses the 


















































function 


0.8phg(x\ .9) + 0.1phg(x; —. 6 ) (25) 

+0.04/i(a:; . 2 , . 01 ,3) + 0.06/ 2 (a:; . 6 , . 02 ). (26) 


Scattering kernel P7 


p(x) = 


where fi and / 2 are as in eq. l |23||24| . The integrand in eq. 
turns out to be analytic in the strip S a with a « 10 2 . The 
scattering kernels Po and P 7 for this phase function were com¬ 
puted at 200 x 200 points with eq. |l9| ), using N — 128 terms 
in each case. A logarithmic heat map of Po is shown in fig. [5] 
and a heat map of P 7 is shown in fig. [ 6 ] The calculation took 
about 10 seconds per scattering kernel on a laptop equipped 
with a dual-core processor running at 1.40 GHz. The relative 
accuracy of each result is about 10 -3 . 


LoglO of scattering kernel P0 



Figure 5: Logarithmic heat map of scattering kernel Po 
for the phase function given by eq. (26). 


4 Conclusion 

A direct numerical integration method using the trapezoidal 
rule has been presented for the evaluation of scattering ker¬ 
nels that arise in plane-parallel radiation transfer equations. 
Its convergence is exponential, and the relation between the 
convergence rate and the domain of analyticity of the phase 
function is explained. The note also presents a closed form of 
the scattering kernel for the Henyey-Greenstein phase func¬ 
tion, in terms of complete elliptic integrals of the second kind. 
The closed form can be used to assess the accuracy of the 
proposed numerical integration scheme. 

Most computational approaches to plane-parallel radiative 
transfer use discretizations based on truncated versions of 
eq. § (Nystrom’s method). However, some problems of this 
form also require the evaluation of intensities from scattered 
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Figure 6: Heat map of scattering kernel Pj for the phase 
function given by eq. (26). 



beams which may be computed from scattering kernels. This 
is where a fast and accurate computational scheme such as 
the one presented here will hopefully be of use. 


5 Appendix: Complete elliptic in¬ 
tegrals 

Legendre’s complete elliptic integrals Ko and Eq of the first 
and second kind are defined as 


Ko(m) 

r /2 ds 

Jo y/l — m sin 2 s 

(27) 


1 r ds 

2 Jo yi-f±fcos s 

(28) 

Eo(m) 

rir/2 - 

— / v 1 — m sin 2 s ds 

Jo 

(29) 


1 r L m m 

2 7o V 2 2 

(30) 


where ra £ C is known as the modulus. Note that usually 
these integrals are expressed in terms of the parameter k 
where m = k 2 , leading to the more common notation 

K(k) = K 0 (k 2 ), E{k ) = E 0 (k 2 ) (31) 

For example, the treatment in [5] is in terms of K , E while 
Mathematica® uses Kq , Eq. These integrals converge for 
k £ C with < 1 and the functions can be continued 
analytically to C minus a branch cut along [1, 00 ). It is known 














(see eq. (19.4.1) in [8]) that 

2mf-K 0 (m) = MEi - K 0 {m). (32) 

dm 1 — m 

Therefore, for a, fd £ C, we can set m — and obtain in 
the case when < 1 

Oi + P 


p 

J yjoi - /3 cos s ds = 2^/a + f3E 0 (34) 

where the principal branch of the square root is used. 

Given real 0 < m < 1, then Eo(m) and Ko(m) may be 
evaluated rapidly using the iterations 
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for n = 0, 1,2,.... Then 


(35) 

(36) 
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lim a n — lim g n — M 

n—too n—>• oo 

where M — AGM(1 : go) is known as Gauss’s arithmetic- 
geometric mean; see [8]. The convergence is quadratic. More¬ 
over, lim n ^oo c n — 0, and the convergence is also quadratic. 
One then obtains the values of Ko and Eo from 

K ° (m> =^i’ E ° {m) = U'-fyA <38) 

Due to the rapid convergence, only a few terms need to 
be evaluated. An alternative fast computation method for 
Eo(M) is given in pQ. 
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